Shock Waves in Weakly Compressed Granular Media 
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We experimentally probe nonlinear wave propagation in weakly compressed granular media, and 
observe a crossover from quasi-linear sound waves at low impact, to shock waves at high impact. 
We show that this crossover grows with the confining pressure Po, whereas the shock wave speed is 
independent of Po — two hallmarks of granular shocks predicted recently [l|- The shocks exhibit 
powerlaw attenuation, which we model with a logarithmic law implying that local dissipation is 
weak. We show that elastic and potential energy balance in the leading part of the shocks. 

PACS numbers: 



Many disordered materials, including granular media 
[l|-|j], foams [5| and emulsions [O, LZf, lose their rigidity 
when their confining pressure Pq is lowered. In almost 
all cases, the resulting unjamming transition goes hand 
in hand with the vanishing of one or both elastic mod- 
uli [8l4l7|. and consequently, nonlinearities must domi- 
nate when such marginal systems are subjected to finite 
stresses [l|, |l8|, |l9[- For example, soft particles exhibit 
nonlinear rheology near jamming [6|,l20|, even when their 
local elastic and viscous interactions are linear [2l|, \2^ , 
and marginally connected spring networks exhibit nonlin- 
ear elasticity near their critical points 23|, |2j| . In these 
two cases, the vanishing of the elastic moduli is a col- 
lective phenomenon, closely connected to the isostatic 
character of the marginal state J8l4l7| . 

Here we experimentally probe a different scenario 
where local nonlinearities near unjamming lead to van- 
ishing elastic moduli and nonlinear excitations: shock 
waves in granular media [25|. Granular media have fric- 
tional interactions, and these prevent granular media to 
reach the isostatic limit. Consequently, there are no col- 
lective mechanisms leading to vanishing elastic moduli or 
nonlinearities |15l. |26I428| . Nevertheless, when granular 




media unjam as the pressure Pq is lowered, the individ- 
ual contacts weaken due to the nonlinear local Hertz con- 
tact law, which states that for elastic spheres, the contact 
forces / scale with deformations S as f ^ S^^^ 29|, |30| . 
As a result, frictional granular media have a vanishing 
linear response at their unjamming point, and their elas- 
tic moduli and sound speed vanish as P^' and Pq' , 



respectively [lO|, |3l|-l34 1 . 



Recent simulations on frictionless Hertzian media by 
Gomez et al. suggest that sound waves give way to 
strongly nonlinear shock waves near the unjamming 
(Po -^ 0) point m. Three crucial questions remain open, 
as the numerical model of Gomez has no static friction, 
no dissipation and is in 2D. First, realistic granular media 
are frictional: do shock waves also arise for non-isostatic, 
frictional systems? Second, friction also leads to dissi- 
pation — do shock waves survive realistic levels of dis- 
sipation? Third, can such shock waves be excited in 3D 
experiments? 



FIG. 1: Schematic side view (a) and top view (b) of our setup. 
Shock waves are excited in the granular medium by a plunger 
(A) which slides through a circular hole (B) and is impacted 
by a heavy mass (C). Pressure sensors and accelerometers 
(big resp. small squares - both enlarged for visibility) are 
connected to steel wires immersed in the granular medium at 
locations xi and X2- 



To answer these questions, we experimentally probe 
sound and shock waves by impacting a weakly com- 
pressed granular medium with a heavy mass, while mea- 
suring the propagation speed and front shape for a wide 
range of impact magnitudes (Fig. 1). We find a novel 
power law attenuation of the shock waves, which we can 
model with a simple model where local dissipation de- 
pends logarithmically on grain forces. Notwithstanding 
this weak dissipation, our shock waves exhibit the three 
main hallmarks of the numerically observed conservative 
shocks: (i) a crossover from linear waves to shock waves 
when the impact pressure P exceeds the confining pres- 
sure Pq; (a) independence of the shock speed on Pq but 
powerlaw scaling with impact strength; (Hi) a balance of 
kinetic and potential energies in the leading edge of the 
shock waves. The ease with which we can excite such 
granular shock waves suggests that they play an impor- 
tant role whenever loose granular media are strongly ex- 
cited [ii-m. 

Setup — Our setup consist of a large metal container 
(45 X 45 X 45 cm) filled with glass beads (diameter 
D = 3.8 — 4.4 mm) and covered by an 2 kg aluminum 
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FIG. 2: (Color online) Typical pressure and velocity signals 
for high impact amplitude, for Po = 0.9 kPa, at xi = 5 cm 
(black) and X2 = 15 cm (red/dashed), (a) Local pressures 
P{t). (b) Local velocities u{t) obtained by integrating the 
acceleration signal. 



top plate. To excite shock waves, we impact a freely slid- 
ing cylindrical piston (A), (diameter 10 cm, length 18 
cm, mass 3.7 kg), which makes contact with the grains 
through a circular hole in the side of the container (B), 
with a mass (C, 1.2 kg) suspended from a pendulum 
(mass 2 kg, length 1.3 m). We detect wave propaga- 
tion throughout the material via pressure sensors and 
accelerometers hurried in the granular medium; the sen- 
sors are attached to steel wires to ensure their correct 
positioning and orientation. We thus probe the local ex- 
cess pressure (our pressure sensors have no DC response) 
and acceleration at distances xi ^-nd X2 from the impact 
zone (Fig. 1). See Supplemental Material at [URL will 
be inserted by publisher] for experimental details. 

Phenomenology — In Fig. 2 we show typical time 
traces of the local pressures and velocities detected at 
locations Xi = 5 cm and X2 = 15 cm, for a strong impact 
and low confining pressure. The waves take the form of 
fronts with a clearly identifiable leading edge where the 
pressures and particle velocities peak, followed by a long 
tail with a complex structure. Here we will focus on the 
nonlinear regime, where the propagation is set by the 
amplitude. We characterize our waves by the peak pres- 
sures, Pp^j and P^2 J ^^d peak particle velocities, ui and 
U2, at xi £md X2- We determine the front speed Vs from 
the time of travel At, where At is the interval between 
P reaching its 50% value at xi a-nd X2- 

Before embarking on a systematic study of the propa- 
gation of these fronts, we briefly discuss the typical scales 
shown in Fig. 2. First, the typical grain velocities are of 
the order of cm/s, the impact speeds of the plunger are of 
the order of m/s whereas Vg is of the order of hundreds of 
m/s. Second, the grain displacements at xi at the time 
of the pressure peak are of the order of 10 ^m, consistent 
with purely elastic deformations (See Supplemental Ma- 
terial at [URL will be inserted by publisher] for details), 
whereas the total motion of the plunger into the sand for 
such impacts is of the order of a mm. 
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FIG. 3: (Color online) Front speed Vs as a function of peak 
pressure in the front, Pm, for a range of confining pressure Po, 
at xi = 5 cm at X2 = 15 cm. Each data point corresponds 
to a number of runs, and we indicated a single typical 5% 
error bar. Fits (solid lines) are functions of the form Vs = 
C(l + P^/Pff''. For P„ > Po, Vs ~ P'/" (dashed line) 
and is independent of Po . 



The picture that emerges is that upon impact, a rapid 
front is formed, which, for the cases when Vs is larger 
than the sound speed, we will call a shock wave. We note 
that we do not observe the breaking up of our fronts in 
several distinct solitons — in the simulations of [l|, the 
stability of the shocks was attributed to the disorder of 
the granular packings. The shock being nonlinear, its 
speed is set by its amplitude, and the leading edge of the 
shock remains leading, as it outruns whatever happens 
in the tail. The interactions in the leading edge are dom- 
inated by elastic, Hertzian interactions. However, long 
after the shock wave has outrun our system, the plunger 
is still slowly penetrating the sand bed, leading to a very 
long tail, where the vast majority of rearrangements and 
dissipative events take place. 

Propagation speed — To probe the nature of the waves 
excited in this system, we varied the pressure Pq at the 
depth of the sensors from 0.9 to 6 kPa and have deter- 
mined the propagation speed Vs for a wide range of im- 
pact strengths. To compare our data to the theoretical 
predictions of [l|, we plot our data in a so-called Hugo- 
niot plot. Due to attenuation (see Fig. 2), P^^ < P^^ 
and as a measure of the impact strength we use their 
geometric mean, Pm '■— \J Pxi^xi- 

Fig. 3 shows Vs as a function of Pm and Pq, and that 
our data is weh fit by the form Vs = C(l -f Pm/Pff'^, 
which is very close to the more involved analytical for- 
mula for shock waves presented in [l| . 

The main three features of our data are (i) For strong 
impacts, Vs becomes independent of Pq, and Vs scales 



consistently as P, 



1/6 



the dependence on Pm is a hall- 
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FIG. 4: (Color online) Scatter plot of the maximum pressures 
Pxi and Px2 at xi = 5 cm and X2 = 15 cm for a range of 
pressures — same data as in Fig. 3. The black curve is a 
fit to our model (Eq. ((5|-((4|, See Supplemental Material at 
[URL will be inserted by publisher] for details. ) 



mark of nonlinear waves, and the exponent 1/6 is the 
one predicted for Hertzian shock waves, (ii) For weak 
impacts, Vs becomes independent of the impact strength 
— the hallmark of linear waves — but increases with 
pressure due to the nonlinear local interaction law. (Hi) 
The crossover from the linear to nonlinear regime is ex- 
pected to arise when P^ ^ Po, because when Pm ^ Po, 
the Hertzian interaction can be linearized, whereas for 
Pm ^ Pq, linearization fails and shocks arise [ij. Indeed 
we find that Pf grows with Pq. 

All these features are in good qualitative and quanti- 
tative agreement with the earlier numerical findings of 
Gomez et. ah the waves we excite for large impacts arc 
indeed shock waves. 

Attenuation — The peak pressure diminishes whilst 
the shock propagates through the material. In Fig. 4 we 
compare the peak pressures P-^^ and P^^ at xi = 5 cm 
and X2 = 15 cm for the same experiments as shown in 
Hugoniot plot. Strikingly, the attenuation varies signif- 
icantly with impact strength: for the weakest impacts, 
Px2 ~ Pxi ' "whilfi foi' larger impacts the relation between 



P^2 and P-^j is consistent with power law scaling: 



For P > 1 : P2 
For P < 1 : P2 



Pf 

A, 



(1) 

(2) 



XI 2 /-f*! the characteristic 
0.77 ±0.05. 



where P := P/P*, Pi^2 := P^ 
pressure P* w 50 Pa and /3 

We will now determine a simple local model that cap- 
tures this remarkable powerlaw attenuation. We ignore 
disorder and imagine the signal to propagate over a linear 
chain of beads, where each bead attenuates the signal by 
a factor (1 — e), where the local attenuation e depends 
on the local pressure P. After some algebra, we find that 
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FIG. 5: (Color online) (a) Attenuation per contact e{P). (b) 
Scatter plot of Pi and P2 for Po = 1-4 kPa at various locations 
Xi and X2, and corresponding fits to our model. 



we can capture the power law relation between P-^^ and 
P^2 when e(P) has the following logarithmic form (See 



Fig. 5a) 



For P>P^: e(P) = e. In 
For P <P^ : e(P) ^ , 



(3) 
(4) 



where Eg '■= ^ln{/3) D/{x2 — Xi) is a material constant 
— See Supplemental Material at [URL will be inserted 
by publisher] for details. 

A surprising consequence of this model is that it pre- 
dicts that the exponent (3 should depend on the distance 
between X2 and xi- To test this, we have performed 
several sets of experiments where we vary this distance. 



Fig. 



5b shows that the log slope relating P^^ and P^j 



indeed increases for larger propagation distance X2 — Xi 
— these trends are captured by our model without addi- 
tional fit parameters. 

Local Elastic Motion and Energy Balance — As Fig. 5a 
illustrates, even for the strongest impacts we can access, 
the attenuation per particle is not more than 10% — it 
is likely smaller, as our model underestimates the num- 
ber of contacts per unit length in a disordered media. 
Hence, we believe that the interactions in the leading 
edge of the shock wave arc predominantly elastic. Addi- 
tional evidence for this comes from a direct comparison 
between the contact forces and particle displacements; 
these are consistent with purely elastic (Hertzian) de- 
formations (See Supplemental Material at [URL will be 
inserted by publisher] for details), showing that sliding 
and rearrangements do not play a dominating role in the 
leading edge of the shock. 

This motivates us to probe the balance between elastic 
and potential energies in our shock waves. For the con- 
servative shocks studied in the numerical simulations [iJ , 
the kinetic and potential energies balance in the shock 
regime, whereas away from the shock regime, the kinetic 

energy tends to zero, but the potential energy saturates 

5/3 
at its lower bound ^ Pq ■ 

To probe this balance in our experimental data, we 

present scatter plots of the peak velocities ui and U2 vs 

the peak pressures P^^ and P;^^ in Fig. 6. The kinetic 
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FIG. 6: (Color online) A comparison of the peak velocities 
and peak pressures at (a) xi = 5 cm and (b) X2 ~ 15 cm 
shows a striking similarity with the numerical data of Gomez 
at. 



energy simply scales ^ uf^ where u„i is the maxiinum 
local velocity, whereas the potential energy for Hertzian 
contacts scales ^ Sm ~ fm , where 6m and fm are the 
local maximuni deformations and forces — note that the 
contact force fm should contain both the DC component 
oc Pq and AC component oc P^-^ ^ . Our data shows a 
strikingly good qualitative agreement with the numeri- 
cal data by Gomez et al. (Fig. 3d of [l[): it is fully 
consistent with a 5/6 power law in the strongly nonlin- 
ear regime [40| , and a pressure dependent deviation from 
this law in the weakly nonlinear regime. Moreover, an 
estimate of the kinetic and potential energies per parti- 
cle based on the data shown in Fig. 6, shows that they 
are of similar order, just as in the simulations [l[. We 
finally note that this balance is equally good at xi and 
X2, even though the energies in X2 are smaller than in xi 
due to attenuation. All this suggest that the attenuation 
does not significantly upset the balance between kinetic 
and potential energies. 

Discussion — By varying the impact strength, we can 
tune our waves from pressure dependent sound waves 
at low impact, to pressure independent nonlinear shock 
waves at higher impact, similar to what was predicted 
for dissipationless Hertzian particles [1[. We find that 
propagating from one grain to the next, a small amount 
of energy is dissipated, which leads to novel powerlaw 
scaling of the attenuation of the shocks amplitude over 
several decades. Is this attenuation caused by scatter- 
ing or dissipation? While we cannot rule out scattering, 
we note that the overall decrease of the pressure and 
velocity profiles shown in Fig. 2 favor a dissipative pic- 
ture; in addition, the width of the leading edge does not 
change very much under propagation. Notwithstanding 
this weak dissipation, the magnitude of displacements in 
the shock, and the balance of kinetic and potential en- 
ergy, show that the physics in the leading edge of the 
shock is dominated by elastic interactions. 

The shock waves we observe here are therefore qualita- 
tively different from the slow granular dcnsification waves 



known as "plowing" (39|, |4l| . Plowing is associated with 
dcnsification through highly dissipative rearrangements, 
and such dcnsification fronts propagate with velocities 
far below the sound speed [39|. Similarly slow and dis- 
sipative events may occur in the tail of our waves, but 
the point is that the leading edge of our shock waves 
propagate faster than the sound speed. They are also 
qualitatively different from the weakly nonlinear waves 
observed under continuous driving |33l.l34|. 

Finally we note that there are two pressure scales that 
set the physics of granular shock waves in experiments — 
the external pressure Pq that sets the crossover from lin- 
ear to shock waves, and the characteristic pressure P* 
above which attenuation sets in. In our experiment, 
Pq ^ P* , but it is conceivable that for more elastic parti- 
cles, or in microgravity, one can reach Pq <C P*, in which 
case, virtually dissipation free granular shock waves could 
be produced. 

We acknowledge discussions with J. Burton, L. Gomez, 
H. Jaeger, X. Jia, and V. Vitclli, and funding from FOM, 
SHELL and NWO. 
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Supplementary Material 

Experimental Details — To detect the propagation of 
waves throughout our granular material, wc employ sev- 
eral Briiel & Kjaer sensors. First, at each location xi, X2 
we use pairs of accelerometers (Deltatron 4508 B002, 1 
V/g and 4508 BOOl, 10 mV/g, resonance at 26 kHz, di- 
mensions 1x1x1.6 cm, weight 4.8 g). Together these cover 
the wide range of accelerations that we encounter (from 
less than 0.1 m/s^ to larger than 7000 m/s^), although 
they saturate for some of the very strong impacts, which 
we conclude to lead to accelerations exceeding 700 g! The 
accelerometers allow us, by integration, to measure the 
local particles velocities and displacements. 

We also use force sensors (Type 8230, 112.5 mv/N at 
Xi and 82301-001, 22.05mV/N at X2, resonance at 75 
kHz, cylindrical shape 19 mm long, circular contact area 
1.33 cm^, weight 30 g). These do not saturate for our ex- 
periments and we focus on the signals from these sensors 
in most of our analysis. The sensors are much stiffer 
(2kN//im) than the typical stiffnesses of the Hertzian 
contacts we encounter, and are immersed in the grains, 
thus acting as pressure sensors. They have no DC re- 
sponse, so that they only are sensitive to the deviation 
from the confining pressure. To convert force to pressure, 
we divide by their circular contact area. 

In our experiments, we vary the pressure Pq at the 
depth of the sensors from 0.9 to 6 kPa by varying the 
load on the topplate. For each Pq, we perform several se- 
ries of experiments, each consisting of hundreds of impact 
experiments. We vary the impact strength (characterized 
by the peak pressure) over four decades: lightly tapping 
the resting pendulum generates the weakest impacts we 
can detect, while swinging the pendulum with full force 
from its maximum height generates the strongest im- 
pacts — only at the lowest confining pressures and at 
the strongest impacts, the piston penetrated the pack- 
ing significantly, limiting the number of experiments in 
this regime. Before each scries of experiments the con- 
tainer was freshly filled and then tapped to stabilize the 
packing, while after each series the container was emp- 
tied and the correct placing of the detectors was verified. 
No systematic dependence on filling procedure or "age" 
of the sample was detected. For each impact, the signals 
from all six sensors were recorded at 250 kHz, and data 
sets ranging from 300 point before to 2000 points after 
impact were stored. To accurately determine the timing 
where the pressures reach their 50% peak magnitude, we 
perform linear fits to the 20% - 80% leading slope of the 
force signals. 

Wave profiles — Our sensors allow us to extract the 
time evolution of the local grain pressure, acceleration, 
velocity and displacement. Fig. ?? shows representative 
examples of these profiles, both for a weak and a strong 
impact, qualitatively illustrating important aspects of 
the phenomenology of the waves we observe. 

These signals allow us to extract typical scales charac- 



terizing the wave profiles. The first striking observation 
is that even for strong impacts, the displacements are 
rather small, and consistent with Hertzian deformations. 
Assuming that approximately 10 particles make contact 
with the sensor, we find that for the string shock wave 
shown in Fig. ??e-h, the peak force at xi is of order 0.6 N 
(and is reached at time t* ^ 0.05 ms). Using Hertz law, 
which states that the contact force F varies with parti- 
cle deformation S a.s F ^ JE*R^/^S^/'^ ^, and taking 
E* = 50 Gpa (typical value for glass), we estimate that 
6 ~ 0.35 fim. The location xi is 12 particle diameters 
w 24 contact zones away, so the cumulative motion x at 
Xi, assuming a fairly flat pressure profile, is predicted 
to be of order 8^m, which is close to the measured dis- 
placement at t* (8/i m, see Fig. ??h). Such estimates are 
also reasonably good for other nonlinear waves, which 
suggests that in the leading edge of the wave, the de- 
formations are predominantly elastic, and grain sliding 
and rearrangements have not arisen yet (these do arise 
in the long tail of the wave [33|, |3^). The elastic nature 
of the contacts also explains why the pressure and veloc- 
ity signals look qualitatively similar — local forces are 
proportional to the pressure, while local deformations S 
of the grains are given by the gradient of x, which, as- 
suming that the pulse moves with constant speed Vs, is 
proportional to the local grain velocity Ug. 

Model for Attenuation — In our simple model for at- 
tenuation, we ignore disorder and thus imagine the signal 
to propagate over a linear chain of beads, where each 
bead attenuates the signal by a factor (1 — e) — see 
Fig. 8a. We take D as our length scale {x := x/D), 
define S := (x2 — Xi)/D, and in Eq. ([5]-??) below drop 
the tildes for simplicity, so that 

P'+i = (1 - eiP')) P' , (5) 

and after taking the continuum limit we obtain: 

To infer e{P) from the observed relation between P2 and 
Pi (Eq. [l][2]) is far from trivial, as P2{Pi) follows from 
integrating Eq. ^ over a finite distance S. 

To determine e{P) we note that solutions of Eq. ([6|), 
P(x), can all be derived from a single mastercurve P{C)- 
P{x) ~ V{x—dx), where the shift dx is adjusted to satisfy 
the initial condition of Eq. © — see Fig. 8a. Hence, we 
obtain the following relation between V, Pi and P2: 



V~\P2)^V-\Pi) + S 



(7) 



Focussing on the power law relation between P2 and 
Pi for P < 1, we need to determine V~^ that satisfies: 
V-'^{P'i) = r-\Pi) + S , which suggest that V'^ should 
be a double log function — so as to turn the exponent /S 
into the additive term S. After some algebra, wc obtain 
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FIG. 7: (Color online) Typical signals for low (a-d) and high (e-h) impact amplitude, for Pq — 0.9 kPa, at Xi = 5 cm (black) 
and X2 ~ 15 cm (red/dashed), (a) Local pressure P{t), and definition of peak pressures Pi, P2, slope (only for W2), and At 
as delay between midpoints of slopes, (b) Local acceleration x{t). (c) Local velocity 11(f) :— L drxij). (d) Local displacement 
X := L dru{T) - the star indicates the peak time of the pressure, (e-h) Same signals, now for much larger impact. Note: the 
data in panel e is identical to the data shown in Fig. 2a in the main text, but here we take a linear scale and the raw contact 
force instead of pressure. 
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FIG. 8: (Color online) (a) Illustration of particle model and 
relation between 7'(C), ^i, P2 and 5. (b) Actual solution for 

P{x). 



that V{C) is a double exponential function: 



p-i(C) 



s 



ln{l3) 



ln{ln{C)) 



P{C) = exp(exp(^C)) , 



(8) 
(9) 



where we stress that V asymptotes to P* , indicating the 
absence of dissipation when P ^ P* (Fig. 8b) . 

This solution P(C) satisfies Eq. IHlwhen £{P) has a log- 
arithmic form, which in dimensional quantities is given 
by Eqs. ([3]) and (gj). 

Potential and Elastic Energies — The peak kinetic en- 
ergy Uk per particle measured at Xj equals the (l/2)mw^, 
where for average grains (r w 2 mm), ra k, p A/Sirr^ w 



8.7 10 ^ kg. The peak potential energy Up per con- 
tact equals L f{x)dx, where 5 is the maximum defor- 
mation and / the contact force, which we assume to 
be Hertzian: f{5) = K6^/^, where K = 4/3 E*r^/^. 
Hence, Up = 2/5 KS^/"^ = 2/5 i^-2/3j5/3_ Assuming 
that 10 particles are in contact with the force sensor of 
area A — 1.33 10~^ m^, and taking into account that the 
total pressure is the sum of Pq and P^. , we can estimate 
the individual contact forces as / = A/10 {Pq + P^t)- 

Now assuming that the peak potential energy of z con- 
tacts equals the peak kinetic energy, we arrive at 



1/2 muj = 2/5 zii:-2/5/5/3 



(10) 



u. - \[^P^-^^^i^' (Po + Px.f" • (11) 

The straight lines in Fig. 6 correspond to z = 5, which is 
a reasonable number for frictional particles in 3D, which, 
in the limit of hard frictional spheres have 4 < 2 < 6. 



